Fixation in a cyclic Lotka-Volterra model 
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We study a cyclic Lotka-Volterra model of N interacting species populating a d-dimensional 
lattice. In the realm of a Kirkwood approximation, a critical number of species N c (d) above which 
the system fixates is determined analytically. We find N c — 5,14,23 in dimensions d — 1,2,3, in 
remarkably good agreement with simulation results in two dimensions. 
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A cyclic variant of the Lotka-Volterra model of in- 
teracting populations, originally introduced by Vito 
Volterra for description of struggle for existence among 
species has then appeared in a number of appar- 

ently unrelated fields ranging from plasma physics || 
to integrable systems Recently, the cyclic Lotka- 

Volterra model (also known as the cyclic voter model) 
has attracted a considerable interest as it was realized 
that introduction of the spatial structure drastically en- 
riches the dynamics |^11|. Namely, if species live on a 
one-dimensional (ID) lattice, a homogeneous initial state 
evolves into a coarsening mosaic of interacting species. 
This heterogeneous spatial structure spontaneously de- 
velops when the number of species is sufficiently small, 
N < N e , where N c = 5 in one dimension P,p|Jll|]. For 
N > N c fixation occurs, i.e. the system approaches a 
frozen state. Little is known in higher dimensions, even 
existence of N c has not yet been established theoretically 
or numerically (in simulations on 2D lattices with N < 10 
species, no sign of fixation has been found and instead a 
reactive steady state has been observed |j|-|ll]]). In this 
work we investigate the cyclic Lotka- Volterra model in 
the framework of a Kirkwood-like approximation. This 
approach predicts a finite N c in all spatial dimensions. 

In the following, we shall use the language of voter 
model |l2] . Consider the cyclic voter model with N pos- 
sible opinions. Each site of a d-dimensional cubic lattice 
is occupied by a voter which has an opinion labeled by a, 
with a — 1, . . . , N. Voters can change their opinions in 
a cyclic manner, a — > a — 1 modulo N, according to the 
opinions of their neighborhood. Specifically, the follow- 
ing sequential dynamics is implemented: (i) we choose 
randomly a site (of opinion a, say) and one of its 2d 
nearest neighbors (of opinion /?); (ii) if /3 = a — 1, then 
the chosen site changes its opinion from a to f3 = a — 1 ; 
(iii) otherwise, opinion does not change. We set the time 
scale so that in unit time each site of the lattice is chosen 
once in average. When N = 2 the cyclic voter model 
is identical to the classic voter model which is solvable 
in arbitrary dimension fl3|| ; therefore in the following we 
assume that N > 3. 

In order to simplify notations, we consider first a ID 
chain. We define p ai ,...,a i {'t) as the probability that a 
randomly chosen segment of i consecutive sites contains 
opinions «i, . . . , a^. For instance, the one-point function 



p a (t) is just the density of opinion a. It obeys 



' Pa-\-l,a Pa,a — 1 Pa — la, 



(1) 



We consider random and uncorrelated initial opinion 
distributions. This implies p a (0) = 1/N, and gener- 
ally Pai,...,ai (0) = 1/N 1 . Symmetry leads to p a , a +i — 

Pa+l,a = Pa-l,a = Pa,a~l, SO Eq. (|j) gives dp a /dt = 

and hence p a (t) = 1/N. Although the dynamics is non- 
conserved, i.e. the densities can change locally, we see 
that for the symmetric initial conditions with equal con- 
centrations the densities are conserved globally. The two- 
point functions obey 



j dpa,a 

dt 



Pct—X,ct,a Pa , a , a 1 T" Pa-\-l,a 

(2) 



r.dp a ,a+l 

^ 7~. — Pq.q+1 Pa—l.a,a+l Pa.a+l,a 

dt 

+ Pa,a+l,a+l + Pa,a+2,a+lj 

which are valid for arbitrary N > 3, and 

r.d<Pa,a+i 

^ ~ Pa — l,a.a+i Pa,a-\-i,a-\-i— 1 



(3) 



dt 



' Pa.a-\-l.a-\-i ~r Pa,a-\-i-\-l,oc-\-i' 



(4) 



Eqs. (D apply lor N > 4, 2 < i < N - 2. Of course, the 
indices are taken modulo N. Finally, for symmetry rea- 
sons p a ,a+N-i = Pa,a-i = Pa.a+i, and more generally 

Pa,a+N — i — Pa,a+i- 

Above equations are exact and normalization can be 
verified. For instance, Y,i<t<N Pa,a+t = Pa = 1/N. 
Eqs. (|l])-(|4|) are the first of an infinite hierarchy of equa- 
tions which is hardly solvable. However, a considerable 
insight can be gained within the two-sites mean-field ap- 
proximation (also called Kirkwood approximation) that 
expresses fc-point functions via one- and two-point func- 
tions fl4|| . For example, the three-point functions read 



Pa\ ,CK2,C*3 



Pa\ .aiPa^ ,as 
Pa 2 



(5) 



This kind of factorization approximation originally de- 
veloped in the realm of equilibrium statistical mechanics 
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has proven to be remarkably successful for a number of 
non-equilibrium processes as well 15|. 

The ansatz of Eq. (||) closures the above rate equations; 
e.g., Eqs. (||) become 



Nr-i 



(ri_i - 2n + r i+ i) , 



(6) 



where rj = p QiQ +j, so for instance r% = p a , a +i is the 
concentration of reactive pairs. Note that the evolution 
rules which define the model are translationally invariant 
in "opinion space" and therefore for translationally in- 
variant initial distributions the two-point correlator p a ,p 
is only a function of (3 — a. Hence Nri is the probabil- 
ity that opinions of any two randomly chosen consecutive 
sites differ by i. The normalization condition thus reads 



<i<N-l ' i 



try requirement, r, 
yields 



rjv- 



the normalization condition 



M-i 



ro 



ro 



2j> 

8=1 

M 



r M = — , N = 2M, 
N 



N = 2M+1. 



(7) 

(8) 



We now turn to the arbitrary dimension d. Making 
use of the compact notations rj, we arrive at the gener- 
alization to the previous rate equations (valid within the 
realm of Kirkwood approximation) 



2d 
2d- I 

2d 
2d - 1 



Nr 



(2d-l)N 
1 



2d 



(2d - 1)N 
Nri [n-i - 2ri + r i+ i] , i 



2r + 2n 
fr - 2r x + r 2 
2,.. 



, (9) 
,M- 1. 



The last equation looks different for even and odd N: 
2d - 1 



Nn (2r M -i - 2r M ) , N = 2M, 
N = 2M+1. 



2d 
2d -I 

Tm = 2d Nri \r M -i - r M ) 



We have two stationary solutions. The first is 

1 



r l = r 2 = ■ ■ ■ = r M = r 



(10) 
(11) 

(12) 



(2d- 1)N 

which together with the normalization condition yields 
2d -2 1 



ro 



(2d-l)N 2 (2d-l)N 

(2d - 2) 
(2d-l)N 2 



for 



1. 



,JV-1. 



(13) 
(14) 



This solution describes the reactive steady state. Note 
that fj oc (d — 1), implying a drastic difference between 



ID and higher dimensional systems. In ID, fj = cor- 
responding to coarsening is feasible, while for d > 1 we 
have fi > implying to a reactive steady state. The 
second stationary solution 



fi =0, 



7^ when i ^ 1 



(15) 



corresponds to fixation; it is possible in arbitrary dimen- 
sion. 

To figure out which of these two solutions actually ap- 
pears in the long time limit let us solve Eqs. (||). To 
accomplish this we first replace variables t and rj (t) by 



2d 



(16) 







1/N. Upon combining with the symme- an d 



Ro(T)=r (t)- 



:, R i (r) = r i (t), (17) 



(2d - 1)AT 

In these variables, Eqs. (|^) acquire a pure diffusion form 
dRi 



—r— - Rj-i - 2R.j + R j+1 . 



(18) 



In these equations index is defined modulo N as previ- 
ously. Equivalently, we may treat Rj(r) as a periodic 
function of j. The initial condition reads 

R j (0) = \Y~ W ^ W ' ^° (mOdA0; (19) 
[773-, otherwise. 

Solving (jU) subject to (|l|) yields 
1 ,00 

R ^ = - (2d-TjN £ ^WSr), (20) 

3—-00 

where /„ denotes the modified Bcsscl function of order n. 
If the variable R\(t) — r\(t) remains positive, the modi- 
fied time variable r behaves similarly to the original time 
variable i; in particular, fj = r,-(£ = 00) = i?;(r = 00). 
The latter quantity is easily found (from the general 
properties of diffusion equation) to be equal to the aver- 
aged initial value. Thus i?j(oo) = ^d-i)^ ' an< ^ t nere " 
fore we recover the reactive steady state of Eq. (|l3]). On 
the other hand, if Ri(t) becomes equal to zero at some 
moment Tf, this will be the end of evolution as r = rt 
would imply t = 00. This case thus corresponds to fixa- 
tion: fi = 0, fj = Ri(rf) > for other i. 

Practically, it is convenient to determine the minimum 
of Ri(t) in the range < r < 00; if the minimum is 
negative, fixation does happen. It turns out that the 
minimum becomes negative for sufficiently large N. This 
allows us to keep only the dominant term from the infi- 
nite sum (|20|), so 



Ri(r) = 



1 



-It 



h(2r) 



N 2 (2d-l)N 



(21) 
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The minimum is reached at r = t, = 0.77256363, and 
i?i(r») becomes negative when N > 4.564293 x (2d - 1). 
Given N should be integer it implies N c = 14 in 2D. 
Would we keep all terms in the sum, we would get a 
little smaller non-integer threshold but still the same 
N c (2) = 14. This assertion can be checked numerically 
with great accuracy if we note that the sum in ( p0| ) can 
be significantly simplified. Indeed, using the well-known 
identity M 



z^(2T)=exp[( 



z + z r 



(22) 



j = -oo 

one can derive 



N-l 



h+N 3 {2r) = l£ C~ p exp[(C p + C p )r], (23) 



] = -co 



p=0 



with C = exp(2ny^l/N). Combining @ and @ we 
arrive at 



1 1 " ^exp[(CP + r p -2)7 

UJ TV 2 (2d-l)iV 2 ^ 



which involves only finite summation. This expression 
has been used to check that indeed Ri(t) remains posi- 
tive only for N < 14 in 2D. One can compute N c (d) in 
arbitrary dimension; for instance N c — 5, 14, 23, 32, 42, 51 
when d = 1, 2, 3, 4, 5, 6, respectively. 
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FIG. 1. Time dependence of the concentration of reactive 
pairs ri(t) in two dimensions. Shown are Monte Carlo simu- 
lation results for N = 12, 13, 14, 15 (top to bottom). 

Thus we have found the critical number of opinions 
N c (d) within the realm of Kirkwood approximation. To 
determine actual N c , numerical simulations have been 
performed. We have considered 2D square lattices (max- 
imum size 2048 x 2048) with periodic boundary condi- 
tions. On lattice of these sizes, fixation has been found 



for N > 14. The concentration of reactive pairs ri(t) 
is drawn on Fig. 1 for N = 12, 13, 14, 15 (the simulation 
data represent an average over 20 different realizations) . 
Fig. 2 plots the concentration of reactive pairs provided 
by the Kirkwood approximation. 




FIG. 2. Numerical integration of Eqs. (9) in two dimen- 
sions. Shown are the concentrations of reactive pairs for 
N = 12, 13, 14, 15 (top to bottom). 

A word of caution is in order. For large N, the concen- 
tration of reactive pairs saturates at a very small value. 
Statistical fluctuations around this value may drive this 
value to zero, which is an absorbing state. This would 
imply an apparent fixation. Even for sufficiently large 
systems a few samples with N = 13 have reached this ab- 
sorbing state. However, the role of fluctuations reduces 
with size, and for linear sizes of order 256 and higher we 
have typically seen a reactive steady state when N = 13. 
In contrast, fixation has always been observed for N = 14 
for linear sizes up to 2048. Strictly speaking, our numeri- 
cal results provide a lower bound for the threshold value: 
N c > 14. However, present data support much stronger 
assertion N c = 14, identical to our theoretical prediction 
based on the Kirkwood approximation. 

To demonstrate the validity of the Kirkwood approx- 
imation it is instructive to apply it to the cyclic voter 
model in ID where a variety of results were already es- 
tablished §,1,0. For N = 3 and d = 1 we solve rate 
equations to find 



ri(t) =r 2 (t) = 



1 



1 



9 l + t/2 

Similarly, for N — 4 and d = 1 we find 
1 1 

n(t)=r 3 (t) 



(24) 



16 l + t/2' 
1 1 1 



vT+1/2 16 l + t/2 



(25) 



Thus in both cases the Kirkwood approximation predicts 
1/t decay of the density of reactive interfaces. The long 
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time behaviors for N — 3 and N = 4 cyclic voter model 
in ID agree with our previous mean-field results for these 
cases O], Compare to exact results O, however, mean- 
field treatments predict faster kinetics; e.g., the density 
of reactive interfaces decays as i -1 / 2 and t~ 2/>3 for N = 3 
and 4, respectively As for the threshold number, 

both rigorous approaches and mean-field treatments give 
the same value N c (l) = 5. This suggests that N c (d) 
given by the Kirkwood approximation might be exact in 
higher dimensions as well. 

In summary, we investigated the cyclic lattice Lotka- 
Volterra model. We argued that for sufficiently large 
number of species, N > N c , fixation occurs. Within 
the framework of Kirkwood approximation, the thresh- 
old value N c (d) has been found analytically in arbitrary 
dimension; for instance, N c — 5,14,23,32,42,51 when 
d = 1,2, 3, 4, 5, 6. In one dimension this prediction is ex- 
act and in two dimensions it agrees with our numerical 
findings for lattices of size up to 2048 x 2048. 
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